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The flux of ultra high energy cosmic rays (UHECRs) at E > 10 eV is believed 
to arise in plasma shock environments in extragalactic sources. In this paper, we 
present a systematic study of particle acceleration by relativistic shocks, in partic- 



ular concerning the dependence on bulk Lorentz factor and the angle between the 
magnetic field and the shock flow. For the first time, simulation results of super- 
and subluminal shocks with boost factors up to T = 1000 are investigated and com- 
pared systematically. While superluminal shocks are shown to be inefficient at the 
highest energies (E > 10 18 ' 5 eV), subluminal shocks may provide particles up to 
10 21 eV, limited only by the Hillas-criterion. For the subluminal case, we find that 
mildly-relativistic shocks, thought to occur in jets of Active Galactic Nuclei (AGN, 
r ~ 10 — 30) yield energy spectra of dN/dE ~ E~ 2 . Highly-relativistic shocks ex- 
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■ pected in Gamma Ray Bursts (GRBs, 100 < F < 1000), on the other hand, have 

spectra as flat as E -1 ' 5 . The model results are compared to the measured flux of cos- 
mic rays (CRs) at the highest energies and it is shown that, while AGN spectra are 
well-suited, GRB spectra are too flat to explain the observed flux. The first evidence 
of a correlation between the cosmic ray flux above 5.7- 10 10 GeV and the distribution 
of AGN by Auger are explained by the model. Neutrino production is expected in 
GRBs, either in mildly or highly relativistic shocks and although these sources are 
excluded as the principle origin of UHECRs, superluminal shocks in particular may 
be observable via neutrino and photon fluxes, rather than as protons. 
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1 Introduction 



The observation of the energy spectrum of ultra high energy cosmic rays (UHE- 
CRs) indicates the presence of an extragalactic component at E > 10 18,5 eV. 
Active Galactic Nuclei (AGN) and Gamma Ray Bursts (GRBs) seem to be 
the two most promising source candidates for the production of ch arged cos- 



mic ra y s (CRs). Wo rk in the late 1970s by a number of authors, e.g. iKrymskii 



fll977h : Es3 fll978al lbh. who based their idea on the original Fermi accelera- 



tion mechanism, first presented by iFermil (119491 . Il954l ). established the basic 



mechanism of particle diffusive acceleration in non-relativistic shocks. In this 
mechanism, individual particles are accelerated in a collisionless magnetised 
plasma by scattering off magnetic irregularities to recross a shock front many 
times. Since then, considerable analytical and numerical investigations have 
been performed, but more questions need to be answered concerning the ac- 
celeration mechanism at relativistic shock speeds. 

In this work, we will present a series of simulation studies varying the shock ve- 
locity and the magnetic field inclination to the shock normal, applying various 
particle scattering media, fixed in the plasma flow frame, aiming to provide 
a more refined determination of the possible spectra which could result. All 
shocks under investigation are taken to be oblique, so that the magnetic field is 
neither aligned with nor strictly perpendicu lar to the shock front norm al. Par- 



allel shock simulation studies were given in lMeli and Quenbyi (l2003al ). Monte 
Carlo calculations will be performed in the relativistic shock environments 
believed to occur in AGN and GRB jets. The resulting very high energy CR 
spectra and the subsequent multiwavelength radiation observation from can- 
didate sources, such as AGN and GRBs, will be mentioned. 



Previous work by iMeli and Quenby fl2003al lbh. on simulation studies for sub- 



luminal and superluminal shocks did not establish the relationship between 
spectral slope, shock inclination angles, shock velocity and particle scattering 
model which will be considered here. These past studies amongst others in- 
vestigated the 'traditional' large angle scattering model and the pitch angle 
diffusion model of particles, but only for the extreme limiting case of 89 < 1 /Y. 
In the present work we will allow for pitch angle scattering to vary between 
1/r < 59 < 10/r to correspond to a variety of scattering wave models. We 
will establish the spectral index dependence on the Y of the shock. The shock 
obliquity will be varied while using a series of high velocity plasma flows, rang- 
ing from T = 10 — 1000. The contribution of the high energy CRs from AGN 
and GRBs to the observed diffuse CR spectrum will be discussed, based upon 
the simulated spectra. 

In Section 2, details of the simulation are discussed. Section 3 presents the 
resulting spectra emphasising the spectral dependence on the relativistic flow 
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Lorentz factor. The calculated spectra are used to estimate a possible contri- 
bution of AGN and GRBs to the observed diffuse spectrum of charged CRs. 
The implications for high energy neutrino and photon emission is discussed in 
Section 4. Finally, in Section 5 the results are summarised. 



1.1 Source Candidates for relativistic shocks 



The energy budget available for UHECRs production in GRBs is approxi- 
mately 2 x 10 44 erg/Mpc 3 /yr, as concluded from the observed gamma-ray 
GRB release rate. This energy budget corresponds to the required energy re- 
lease rate of > 10 18 5 eV of UHECRs and is based on an assumed star formation 
rate (SFR) history (Vietri, 1995 ; Waxman . 2000l ) with about 1 burst/Gpc 3 /yr 
at z — 0. The typical luminosity of AGN, L ~ 10 42 — 10 47 erg/s is also suffi- 
cient to explain the UHECR flux under the assumption that AGN follow the 
SFR. For example, iMoran et al.l (120011 ) find an average luminosity from an 
estimated 95 X-ray active AGN within 60 Mpc of 4.8 x 10 44 erg/s, this dis- 
tance being the cosmic ray absorption horizon at about 2 x 10 20 eV. Within 
the regions of local space accessible to cosmic ray diffusion, the energy supply 
over a Hubble time is 6.6 x 10~ 16 erg/cm 3 . By comparison the GRB supply 
is 6.7 x 10~ 20 erg/cm 3 . With GRBs as the most luminous transient objects 
in the sky and AGN as the most luminous permanent ones, these two source 
classes are th e best c andid ates for the acceleration of UHECRs, following the 
arguments of iHillasi (119841 ) . 



The observation of electron synchrotron radiation in the radio regime indicates 
that mildly-relati vistic shocks of boost facto r s r ^ [10 — 30 are pres ent in 



the jets of AGN (IBiermann and Strittmatterl . 119871 ; iFalcke et al.l . Il995l ). The 



photon spectrum of AGN is broadband, ranging from radio up to TeV emission 
in the case of optically thin sources. Assuming that hadrons are accelerated 
along with the electrons in the jet, AGN are good candidates to be responsible 
for at least a significant fraction of the extragalactic component of the CR flux. 



The observation of the highly variable, prompt GRB photon spectra at soft 
photon energies (E 1 > 100 keV) indicate the acceleration of el ectrons in highly 
relativ istic shocks of boost factors around 100 < T < 1000, see lHalzen and Hooper 
(120021 ) and references therein. However, mildly relativistic internal shocks may 
also occur within the GRB plasma flow. Pro t ons ar e believed to be accelerated 
up to ~ 10 21 eV as discussed by IWaxmanl (120001 ). The total release of elec- 
tromagnetic energy by GRBs has led to the suggestion that t he CR spect rum 
above the ankle (E > 10 18 5 eV) can be explained by GRBs (jVietril . Il9951 ). 



3 



1.2 Maximum energy 



A basic physical limitation to the maximum energy of the accelerated p articles 
i s the s ize of the acceleration region and the magnetic field present. iParker 
fll958al lblfl Il966h first d i scusse d this CR energy limitation within the solar 
system. Later on iHillasI (119841 ) extended the CR energy limit argument to 
a number of astrophysical sources. Particles must escape once the gyration 
radius exceeds the source radius. The maximum energy is then given as 



E, 



18 



A • Z ■ B^g ■ Lk P c ■ 



(1) 



Here, -E^fax := -E max /(10 18 eV) is the maximum energy that can be achiev ed, 
A = V s /c and is equal to 1 for the oblique shock conditions (jJokipiil . 119871 ). Z 
is the charge of the accelerated particle in units of the charge of the electron, e. 
Furthermore, B^g '■= B/(l fiG) is the magnetic field of the acceleration region 
in units of 1 /iG and Lk pc '■= L/ (l kpc) is the size of the acceleration region in 
units of 1 kpc. As discussed by IHillasI (119841 ). AGN cores, radio galaxy lobes, 
and hot spots, are promising candidates for the acceleration of the highest 
energy events. A second criterion of Hillas that the proton synchrotron loss 
time should not be less than the acceleration time, is easily met by proton 
shock acceleration if the shock is relativistic since the acceleration time is 



_ 3 x 10 3 XE 18 
i~acc — ^ year , 
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where the scattering mean free path is X times the Larmor radius and the 
synchrotron loss time is 



T, 



synch 



1.4 x 10 14 



year . 



(3) 



As mentioned before, boost factors in AGN are typically assumed to be around 
T ~ 10, although in special cases, the boost factor can reach up to T ~ 30 — 40. 
It is, however, possibl e that the shock itself happens w ithin a general rel- 
ativistic environment ( Biermann and Strittmatterl . Il987l ) and so the shock 
is not necessarily relativistic. Here, we assume relativistic movement of the 
shock front. The magnetic field can be up to B ~ 10 -3 G in the jet at a 
radius of r ~ 1 kpc and the field typically decreases inversely with the ra- 
dius. This condition all ows for particle acceleration up to t he highest energies, 
i.e. i^max^ ~ 10 21 e V ( iBiermann and Strittmatterl . 119871 ). Radio galaxy hot 
spots allow for acceleration of particles up to 10 21 eV and energies of 10 19 eV 
can be produced in radio galaxy lobes. Other suggested sources are not able 
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to produce particles of sufficient energy because the field strengths available 
are insufficient to contain the particles. 



GRBs may also accelerate protons to 5 x 10 20 eV. This is because fields up to 
YB = 1 14 G can be produced at an accretion torus at about 6 x 10 6 cm radius 
in massive sta r collapse before quantum me chanical dissipation, limits the field 
amplification (ILerche and Schramm! . 119771 ) and the subsequent 1/r field fall- 
off yields an E max 
synchrotron losses, 



sec 



Vietri 
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1995 
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Waxman ( 



particles during the prompt emission phase in highly-relativistic shocks (r = 
100 — 1000) is discussed in the following. Acceleration in external shocks during 
the slowing afterglow phase is also possible so that the spectra may resemble 
the AGN spectra that we calculate and present in Section 3. 



2 The physical concept and the Monte Carlo simulations 



Begelman and Kirk fll990h have claimed that most upstream field configu- 



rations at high shock boost factors T appear superluminal. In superluminal 
shocks the particles are accelerated in a shock drift because there is no trans- 



formation into a de Hoffmann- Teller (HT) frame (jde Hoffmann and Teller 



1950|) , where E = 0. To transform from the normal shock frame (NSH) to the 
HT frame, we need to boost by a speed Vht along the shock frame, where 
Vht = Vnsh ■ tan ip (if) is the angle between the magnetic field and the shock 
normal). Due to physical causality, this transformation is only possible if Vht 
is less or equal to the speed of light. Thus, when Vnsh = c, the limit is 
tanip = 1. When tan-^ < 1 the subluminal shock transformation case is ap- 
plied. For all other cases, where a HT frame cannot be found due to a very high 
inclination in combination with a high shock velocity, the superluminal shock 
condition applies. While in the subluminal case particle transmission at the 
shock can be decided in the HT frame employing conservation of the first adi- 
abatic invariant, in the superluminal case computations are followed entirely 
in the fluid rest frames with reference to the shock frame simply employed to 
check whether upstream or downstream shock conditions apply. 

Superluminal shock acceleration may be treated as a shock drift mechanism in 
the shock frame and is best visualised when the shock is nearly perpendicular. 
As viewed in the shock rest frame, the particle is moving in a steep magnetic 
field gradient perpendicular to the shock surface. The plasma motion at an 
angle to the magnetic field creates an electric field given by E = —ux B/c. 
In this mechanism, scattering is unimportant and particles can be accelerated 
by just one shock encounter. However, due to the high field inclination, a 
particle gyrating around the field lines crosses the shock more then once in 
one shock encounter and thus obtains significant additional energy. The drift 
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is in a direction perpendicular to the magnetic field gradient (shock normal) 
and the magnetic field, according to 



- p-c-u.Bx.VB 



with B = \B\. The drift direction is such as to cause the particle's energy to 
increase. 



Meli and Quenby fl2003al lbh showed that a transformation from an initially 



isotropic rest frame distribution to an accelerated flow frame leads to a co- 
moving relativistic plasma frame field distribution lying close to the flow vec- 
tor. This conditio n allows for a range o f sub luminal situations when viewed in 



the shock frame. iBegelman and Kirk! (119901 1 pointed out, on the other hand 



that in the blast wave frame the turbulence can be isotropic and many shock 
stationary frame configurations can be superluminal. In general, flow into and 
out of the shock discontinuity is not along the shock normal, but a transfor- 
mation is possible in t o the NSH frame to render the flows along the normal 
( Begelman and Kirk . 1990l ). We assume that such a transformation has al- 



ready been made. 



Vietri et al.l (120031 1 raise the question as to the correct reference frame in 
which jet acceleration should be viewed. To become CRs, the protons escape 
sideways out of the narrow beamed jet and are seen in the extragalactic rest 
frame. Escape upstream into the ambient medium is very unlikely. Strictly 
speaking, we need to employ a transformation conserving t he distr i bution 



function / depending on the momentum p, so f(p) = f( p >) (iFormanl (Il970h 



demonstrates this invariance which is related to Liouville's theorem). Here, 
P 2 f{p) = dJ(E(p),n)/dE, where f(p) is in the direction of vector n and the 
energy E corresponds to the momentum vector p. The particle flux dJ / dE is 
typically in units of cm" 2 s" 1 sr" 1 GeV" 1 . Prime and unprimed denote rest 
frame and observer's frame. For our calculations we choose the total energy 
normalisation to be done in the shock frame. While the transformation affects 
the maximum energy obtained in the rest frame to some extent, the fact 
that escape from the beamed jet is mainly by motion perpendicular to the 
flow, means that we are chiefly involved in transforming a momentum vector 
perpendicular to the relative velocity of the reference frames, where there is 
no Lorentz correction. For X-ray production in shocks the relevant results are 
in the downstream frame. 

The purpose of the Monte Carlo simulations is to find a solution to the particle 
transport equation for highly-relativistic flow velocities. The appropriate time 
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independent Boltzmann equation is given by the following 

dx dt 



T(V + vfi)^- = %\ c , (5) 



where a steady state is assumed at the shock rest frame, and V is the fluid 
velocity, v the velocity of the particle, T the Lorentz factor of the fluid frame, 
fi = cos6 the cosine of the particle's pitch angle 9 and df/dt\ c the collision 
operator. 

Particle scattering by magnetic irregularities fixed in the plasma frame will be 
assumed. The Alfven velocity, Va, for AGN jets with a maximum 10~ 3 G field 
and a particle density of 10~ 2 cm -3 is 2 x 10 9 cm/s. For GRB jets with a 10 G 
field and a particle density of 10 6 cm" 3 , it is of similar magnitude. With these 
wave velocities, the waves will appear almost stationary to relativistic particles 
and second order Fermi acceleration, with a fractional energy gain (Va/c) 2 , 
is small compared with shock acceleration where a single encounter cause 
energy enhancement by a T 2 factor. Va — > is the theoretical relativistic 

limit to wave speed and to shock speed in the fluid frame while numerical jet 
simul ation suggests large scale, knot features travel at about c/10 (Ivan Putter! 
20051 ). Thus, we seem justified in neglecting fluid frame acceleration beyond 
the region of trajectory intersection with the shock surface. This is especially 
true downstream where the second order Fermi fractional gain per collision is 
limited to 1/12. 



The scattering operator will be treated vi a a pitch angle sca t tering a pproach, 
rather than as previously in the work of Meli and Quenby ( 2003a , bl) where 
large angle scattering was considered. In standard kinetic theory the spatial 
diffusion c oefficien t s k \\ and kj_, are related to the formula kj_ = k\\ ■ (1 + 

\2\ 



Jokipiil (119871 ). In the well known Bohm Limit, A/rj 



1, but in- 
terplanetary particle propagation studies and gyroresonance theory suggest A 
to be a number of particle gyroradii rj, that is A > 10r/. Hence, in the shock 
normal, or in the x direction, the diffusion coefficient is given by k\\ = Xv/3 
where k = n\\cos 2 ip, since we assume that k\\ » k±. A guiding centre ap- 
proximation is therefore used to follow propagation along field lines. Because 
relativistic shocks generate strong small-scale tu rbulent magnetic field down- 
stream by the relativistic two stream instability (jMedvedevl . Il999l ) . we assume 
power in the scattering increases, as does the field strength and for simplicity 
we assume these quantities change to keep their ratio constant. Therefore, 

^down A U p. 



Gallant et al.l (119991 ) have demonstrated analytically that particles entering 
the upstream region in a direction nearly normal to the shock can only experi- 
ence small-angle scattering (pitch angle diffusion), 58 < 1/T with 56 measured 
in the upstream fluid frame for scattering in a uniform field or a randomly 
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oriented set of uniform field cells. The condition arises because particles at- 
tempting to penetrate upstream from the shock are swept back into the shock 
before they can scatter far from it. If we consider the ratio of initial energy to 
final energy, measured in the upstream or 'primed' frame, after up to down to 
up transmission in an inclined shock, we have 

E'/El = T 2 (l + - Pr^u) (6) 



where j3 r is the relative velocity of the frames as a fraction of c, fj,'_^ d is the 
cosine of the particle velocity angle to the normal direction in the upstream 
frame at the moment of going downstream and fi*-> u is cosine of the parti- 
cle velocity angle to the normal in the downstream frame at the moment of 
going upstream. For down to up crossing, —fi*^ u > B\jB\ where B\jB\ is 
the inverse compression ratio. Up to down, the pitch angle s cattering con- 



straint is « -1 + 1/r 2 . As pointed out bv iBarind fll999f ). E' f /E[ « 2. 
On the first up to down to up cycle with injection upstream and directed 
toward s the shock, E'jEj pa T 2 , corresponding to the large angle scattering 



case of iQuenby and Lieul (119891 ) . Since this constraint is largely dependent on 
the kinematic competition between upstream particle flow and the relativistic 
approach of the shock front, it is not critically dependent on the exact mag 



nitude of the pitch angle scattering. Moreover, as shown by IQuenby and Meli 



(120051 ). blobs of high field scattering centres which could cause even large 
angle scattering, are allowed if the Larmor radius in the blob, r 9 .&, satisfies 
(A/2r 2 ) > r g ^ where the parallel mean free path A = Xr g ^mbient and X ~ 10. 
Such fields would allow for substantial scattering before particles are swept 
back to the shock. This criterion, derived for parallel shocks, effectively applies 
in the oblique, subluminal case since in the upstream frame, the field direction 
is near the normal direction. To obtain a deflection ~ lOr -1 , the constraint 
is relaxed to b/B >m 2T as the ratio of 'blob' to ambient field, where b is 
the perpendicular perturbation to the mean field, B. Hence, a relatively few 
field blobs, perhaps originating in strong instability in a hypernova collapse 
and with field strengths up to a factor 1000 stronger than the ambient, would 
thus be required to allow for pitch angle scattering greater than 58 = 1/r. 
For practical purposes, remembering it is the downstream scattering that is 
relevant to particle loss, it seems reasonable to use larger values of 58, within 
the pitch angle diffusion model. 

A simple representation for the effect of the turbulence which relates to previ- 
ous work is to suppose the particle scatters 58 = N/T every A = 10r; where r/ 
is the Larmor radius, in the plane of gyration. A transverse field pertu r batio n 
changes the pitch angle in a quasi- linear theory ( Kennel and Petscheckl . Il966l ). 
by 

58 = uj^-5t (7) 
B 
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in a time 5t due to a perpendicular perturbation, b, to the mean field, B, 
with cyclotron angular frequency, oj = eBj^m c. The particle moves in near 
gyroresonance with the wave in b. A pitch angle diffusion coefficient can then 
be derived 

01 V\\ B 



where P(k) = P k s is the power spectral density of b at gyroresonance wave 
number k = co/vn. A particle then diffuses in pitch by a finite amount, 59 
during St given by 

59 2 = 2D e 5t . (9) 



It is waves, of wave number k, fulfilling the gyroresonance condition, k = oj/v\\ 
that cause the scattering of particles satisfying uri = v±. We choose 59 to lie 
between 1/T and 10/T, so on average a particle scatters 5/T after a time 
10y/3ri/c so that 



2 25 2u 2 y i0V3v± 
56 = — = ——P {k) — . 10 

l z v\\B z vf, ecu 



To choose 59 independent of particle 7 that is of u, we require s=-l for the 
spectral slope. Then we obtain a power spectrum relative to the mean field 
power 

^ - -^k-' (11) 

The total fractional power in the turbulence if resonating waves are present to 
scatter particles of between 7 = 300 and 7 = 10 12 is 1.4/T 2 . Hence, the cho- 
sen pitch angle scattering model corresponds to a weak turbulence situation. 
Because quasi-linear theory for wave particle interactions is known to be an 
inexact approximation, the power spectrum we have presented must also be an 
approximation to the scattering model we employ. The choice of a fixed factor, 
10, in the relation between A and r\ acknowledges this inexactness. The spec- 
trum is simply presented to p rovide a link with the work of others, especially 



Niemec and Ostrowskil (120051 ) who realise the field fluctuations employing a 
specified wave spectrum. In the present investigation, we allow the particles 
with pitch angle chosen at random to lie in the range of 1/T < 59 < 10/T. In 
our past work mentioned in Section 1, we presented results for the large angle 
case and the aforementioned extreme limiting case of 59 < 1/T. 
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Our past studies of pitch angle scattering iMeli and Quenby fl2003al lbh sug- 
gested that as the magnetic field inclination angle to the shock normal de- 
creased, the spectra become smoother. Both a pronounced, plateau-like struc- 
ture and increasing flatness developed for the highest values of the shock boost 
factor T. Additionally, for all inclination values used in the simulati o ns, fo r 
r = 10 — 30, the sp ectral form remains smooth. lEUison and Double! (120041 ). 
Stecker et al.l (120071 ) and references therein, have shown similar trends. In this 
more comprehensive study, the previous claims will be more thoroughly inves- 
tigated. 

Standard theory poses the conservation of the first adiabatic invariant in the 
HT frame in order to determine reflection or transmission of the particles. 
Since in this frame the allowed and forbidden angles for trans mission depend 
only on the input pitch and phase, not on rigidity, the results of lHudsonl (119651 ) 
apply in our model. For an isotropic flux, the transmission coefficient £ is 
simply given by the particle flux conservation between an upstream magnetic 
flux tube area and the corresponding downstream flux tube area to which it 
connects 



FBi 



(12) 



where F is th e distribution function ( jParkerl . Il965l ). Trajectory integration, 
Hudson! (119651 ) giving the phase dependence of the probability of transmis- 
sion as a function of phase and pitch angle showed that the reflection percent 
plotted against pitch angle never varied more than 20% from the mean value 
an d these r esults were consistent with the flux conservation prediction based 
on IParker! (Il965h and the adiabatic invariant conservation. In the relativistic 
shock situation anisotropy renders the input to th e shock from upstream, very 
anisotropic in pitch angle, but as was discussed in iMelil (120031 ). it is an accept- 
able approximation to randomise phase before transforming to the HT frame 
and then to use the adiabatic invariant to decide on reflection/transmission, 
because of the Hudson result. 



For further details on the simulations, see appendix A. 



3 Results 



The physical concepts and analytical approximations previously mentioned 
are used to perform Monte Carlo simulations for relativistic superluminal and 
subluminal shocks. Computed particle spectra will be presented with special 
focus on the relation between a given astrophysical source shock boost factor 
T and the resulting spectral features. Of particular relevance is Section 3.2 
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devoted to the calculation of the contribution of these sources to the observed 
diffuse CR spectrum. 



3.1 CR shock acceleration spectra 



3.1.1 Superluminal shock spectra 

Initially, we present simulations for relativistic superluminal shocks as de- 
scribed in the previous section and the Appendix. Unlike our previous work 
which was confined to large angle scattering, pitch angle scatter is employed. 
The guiding centre approximation is used except within one complete particle 
helical cycle from the shock. Since a transformation into the HT frame is not 
possible as described in Section 2, the particles are followed in the appropri- 
ate fluid rest and SH frames, simulating the physical picture of the shock drift 
mechanism. We follow the helix trajectory of the particle until it intersects 
with the shock front, applying a pitch angle scatter [69 < 1O/T,0 £ (0, 2ir)] 
right up to the shock interface. Simulation runs performed showed that the 
results were basically independent of x/i, magnetic field to shock normal angle. 
We therefore show here a simulation run with an arbitrarily chosen inclina- 
tion angle, ip = 76°, as an example, employing a range of boost factors. The 
resulting spectra are presented in the shock frame on the downstream side as 
we want to ensure comparability with following calculations. 
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Fig. 1. Superluminal, relativistic spectra at if) = 76°. Boost factors are varied be- 
tween r = 10, 100, 300, 500, 1000. Spectra for different inclination angles ip are com- 
parable. 



The resulting particle spectra for T = 10, 100, 300, 500 and 1000 are displayed 
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in Fig. 1. These plots indicate that for mildly- relativistic shocks (T = 10), the 
acceleration is efficient up to energies E p < 10 2 GeV. More highly- relativistic 
shocks (r > 100) can produce particles up to E p < 10 5 GeV. The upper 
limit of significant acceleration is ~ T 2 , corresponding to only one complete 
particle crossing cycle with the majority of particles either failing to return a 
second time or only returning at angles close to the normal. In the region of 
efficient acceleration, the spectra approximately follow power-laws with spec- 
tral indices lying between ~ 2.0 — 2.3. In contrast, for the case of large angle 
scattering, previously studied by lMeli and Quenbyl (j2003bl ). the spectra could 
not be described by power-laws, but exhibited a concave shape terminating in 
a steep energy cut-off. 



We conclude that superluminal, relativistic shocks are not efficient accelerators 
for very high energy particles and are unlikely to contribute to observable 
effects, discussed in more detail in Section 3.2. These conclusions concur with 



the work of iNiemec and Ostrowskil (120071 ). While superluminal shocks cannot 
contribute to the observed spectrum of charged UHECRs, a contribution to the 
neutrino- and TeV-photon background arising from proton-photon or proton- 
proton interactions, is still possible as discussed in Section 4. 



3.1.2 Subluminal shock spectra 

Particle spectra produced in relativistic subluminal shocks have been calcu- 
lated for three different inclination angles in the shock frame, ip = 23°, 33° and 
43° which we chose as representing the possible range of subluminal shock 
angles likely in astrophysical sources. In appendix B, simulated particle spec- 
tra are presented for the entire energy range of particle energy considered, 
E = 10 2 - 10 12 GeV, for a range of boost factors, T = 10 to T = 1000 and for 
all three angles. The chief result is that the spectra appear as smooth power- 
laws for mildly- relativistic shocks, T < 30. As the boost factor increases, the 
spectra appear with bumps (T > 100): The plateau-like parts of the spectra at 
higher energies and higher boost factors are caused by particles continuing to 
undergo significant acceleration in a second cycle. The lower energy part of the 
spectrum is dominated by particles undergoing one acceleration shock crossing 
while the second bump in the spectrum represents particles experiencing two 
acceleration shock crossings. In this first complete shock cycle crossing from 
upstream to downstream to upstream, the energy gain is a factor T 2 , operating 
an injection energy already ~ T in magnitude. Subsequent crossings become 
smoother since the energy gain is expected to be limited to ~ 2 and there is 
a statistical smoothing of the energy gains. At the highest boost factors we 
investigate, this smoothing regime is not reached. 

In order to get a representative picture, we average the particle spectra over 
the three angles at a particular T, in all following calculations. This should 
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give a more realistic view of the diffuse particle flux from extragalactic sources, 
as it is expected that a range of angles is likely to occur in the class of AGN 
and GRB shocks. We concentrate on the highest energies because observation 
of particle-induced air showers at energies between 10 9 ' 5 GeV and 10 10,5 GeV 
indicates an extragalactic origin of the charged CRs, distinct from a dominant, 
galaxy produced component at lower energies. A power-law fit is made to 
the simulated spectra between 10 9,5 GeV and 10 10 ' 5 GeV. Figure 2 shows 
the averaged, simulated spectra between 10 8 ' 5 GeV and 10 11 GeV. At even 
higher energies, the spectrum is altered by the absorption of protons due to 
interactions with the cosmic microwave background. While the normalisation 
of the spectra is arbitrary since dependent on the number of injected particles 
in the Monte-Carlo simulation, the spectral index can be compared to what 
is observed in CRs in the same energy range. Table 1 shows the variation of 
the spectral index with the boost factor. While mildly-relativistic shocks show 
indices around a p ~ 2, highly-relativistic shocks with T > 100 have flatter 
spectra between 0.7 < a p < 1.5. Particle spectra emitted during the prompt 
phase of GRBs (T > 100) will therefore appear much flatter than AGN particle 
spectra (T ~ 10). This has important implications for the interpretation of the 
origin of the UHECR spectrum as will be discussed in a following subsection. 



Stecker et al.l (120071 ) investigating parallel shocks up to T = 30 found increas- 
ing spectral structure and decreasing slope as T increased (i? -1 - 26 at T = 30), 
trends we have shown to extend to far higher r factors and for a more g enera l 



class of subluminal shock inclination angles. iBednarz and Ostrowski fll998h 



used pitch angle scattering and varying cross-field diffusion and found that at 
low T, steep spectra occurred at large inclination angles but all values of these 
parameters seemed to pro duce spectral slopes of -2.2 at T = 243. In contrast, 



Meli and Quenbyl (]2003al ) found spectra flatter than E 1 for parallel shocks 
as T — > 1000. Later work, h owever, with wave spectra P(k) ~ A; _1_>1 ' 5 by 



Niemec and Ostrowski! (120051 ) found spectra flatter than E 1 with noticeable 



structure spectra in inclined, subluminal shocks at upstream velocities of 0.5c 



Their trajectory integrations took into account cross field diffusion. iBaring 



( 120041 ) however, cites previous work with no cross field diffusion which finds 
significant relativistic, inclined shock acceleration to be limited to inclination 
angles < 25°. It is not clear whether the very steep spectra quoted at higher 
inclination angles are due to the steep slope at the edge of the first plateau 
we mention above, or whether there is a significant difference from our mod- 
elling. A high perpendicular diffusion coefficient, high inclination shock situa- 
tion might be expected to approach a high scattering, parallel shock regime, 
but this possibility does not seem to reconcile the various conflicting results 
just discussed. There is an agreement that at very high inclinations significant 
acceleration above that due to a single shock cycle is ruled out. Also, there 
seems to be a developing consensus tha t high T, subluminal shocks result in 



flatter spectra than E 2 . Dingus! (jl995l ) provides gamma-ray burst evidence 



for relativistic electron spectra with relatively flat slopes with exponents at 
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Fig. 2. Subluminal spectra averaged over the three angles for different T: 
r = 10,20,30 is displayed in the first row, in the middle, T = 100,300,500 is 
shown and T = 700, 900, 1000 is the bottom row. The black crosses in each graph 
represent the simulation result. The straight line shows the single power-law for 
comparison. 



least as low as 



-2. 



Results indicating a variety of possible slopes are also consistent with radio 
data on the electron spectra injected at terminal hotspots in the lobes of 
powerful FR -II radio galax i es wh ere n o single, universal pow er-law is found, 
as shown by iRudnick et al.l (119941 ) and iMachalski et al.l (120071 ) among others. 
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r 


10 


20 


30 


100 


300 


500 


700 


900 


1000 


a p 


2.1 


2.0 


2.1 


1.2 


1.5 


1.5 


1.5 


0.7 


1.4 



Table 1 

Spectral indices for a single power-law comparison for subluminal shocks. The spec- 
tral fits were made between 10 95 GeV and 10 10 ' 5 GeV in order to be comparable to 
the observed CR flux at the same energies. The uncertainty from the fit is less than 
10% if we assume an accuracy of the simulation is better than A \og(dN / dE p ) ~ 0.5, 
which is a conservative estimate. 

3.2 Diffuse CR spectra from GRBs and AGN 

The source spectra derived previously can be translated into an expected dif- 
fuse proton flux from astrophysical sources by folding the spectra with the 
spatial distribution of the sources. In this section, AGN and GRBs are used as 
potential candidates because these are the sources with the highest observed 
output in relativistic electrons. Since the particle spectra are strongly depen- 
dent on T, it is important to discuss which spectra to use for these two source 
classes. Spectral choice is investigated in the next subsection before the actual 
calculation of the diffuse spectra is shown. In the last subsection, the results 
of our calculations are compared to CR data. 



3.2.1 AGN and GRBs - intrinsic spectra 

Spectral fits limited to the energy range 10 9,5 GeV to 10 10,5 GeV of extra- 
galactic origin are employed for both AGN and GRB sources, but the boost 
factors applicable differ between these source classes. 

The boost factor deduced from electron synchrotron observation can vary 
significantly in the case of GRBs. While the majority of sources are estimated 
to have boost factors around V « 300, more moderate values down to V = 100 
or more extreme values up to T = 1000 are believed to occur. However, the 
exact distribution of Gamma Ray Burst T factors cannot be determined. In 
many cases, only upper limits can be given. In addition, there may be hidden 
bursts not observed with GRB satellite experiments. It is therefore not useful 
to model a detailed distribution of boost factors for GRBs, while the simulation 
results connecting the spectral index of the spectrum with a boost factor are 
subject to uncertainty, as is implied by the absence of a monotonic trend in 
Table 1. 

It was shown that for boost factors of T > 100, the source spectra lie between 
Ep' 1 ' 5 and E p ~ ' 7 . A conservative estimate for the flattened GRB spectra will 
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be adopted using 

d§GRB 



dE n 



oc E„ 



•1.5 



(13) 



The situation is simpler for AGN, as the shock's boost factors vary only up 
to between T = 10 and T = 30 (IBiermann and Strittmatterl . 119871 ). Here, 
all computed spectra cluster around a value of E P ~ 2A . Therefore, the AGN 
spectrum will be taken as 



d&AGN £ -2.1 

dEn 



(14) 



3.2.2 From CR shock acceleration spectra to a diffuse spectrum 

The diffuse spectrum as measured at Earth depends on several factors: 



Single source spectra at the source d§/dE p . The spectral behaviour was 
already discussed in the previous subsection. To account for particle propa- 
gation, adiabatic energy losses need to be considered as E p (z) = E p ■ (1 + z). 
Here, E p (z) is the energy as observed at a source at redshift z and E p is the 
corresponding energy observed at Earth. Diffusive propagation in the mag- 
netic field between clusters is assumed to involve only small angle scattering 
with preservation of spectral shape. Anisotropy in the source distribution is 
neglected. In addition, we consider pure proton spectra so that spallation 
effects are not present. 

Source evolution g(z): It is assumed that both AGN and GRBs follow the 
SFR to deter mine the numb e r dens ity evolution with comovi ng volume, see 
for ex ample, lHasinger et al.l (120051 ) in the case of AGN and iPugliese et al. 
( 120001 ) in the case of GRBs. A large s ample of radio quiet A GN selected at 
X-ray wavelengths was investigated by lHasinger et al.l (120051 ). The comoving 
density dn/dV(z) is given as 



dn 
dV 



(z) oc < 



1 + z) m for z < zi 

[l + zi) m for z 1 < z < z 2 

[1 + z x ) m ■ 10 fc -( 2 " 22 ) for z > z 2 , 



(15) 



with the parameters m = 5.0, z\ = 1.7, z 2 = 2.7 and k = —0.43. The total 
redshift evolution g{z) further includes multiplying the comoving volume 
dV/dz with a factor l/(Aird 2 L ) to account for the decrease of the flux L 
with the distance di, neglecting a possible travel limitation to the distance 



1(3 



reached by significant magnetic scattering. Therefore, 



dn 
~dV { 



dV 

dz 



(And 



(16) 



For simplicity, this model is used for both AGN and GRBs. Although de- 
viations between the SFR scenarios of AGN and GRBs are expected, the 
approximation that both follow the distribution of radio quiet X-ray AGN 
is reasonable (IHasinger et al.1 . 120051 ): the deviations being expected to be 



negligible with respect to general uncertainties arising from assumptions 
about the acceleration region. 

Absorption of protons at the highest energies: Protons at E p > 5 • 10 19 eV 
are absorbed due to interactions with the cosmic microwave background a s 
was recently confirmed by the Auger experiment (jYamamoto et all 120071 ). 
Therefore, the diffuse spectrum resulting from the propagation of a single 
source spectra is modified by a further factor, exp[— E p /(5 ■ 10 19 eV)] to 
account for this effect. 

The normalisation of the diffuse spectrum: Because the calculated particle 
spectra are given in arbitrary units, normalisation of the overall spectrum 
as measured at Earth is achieved using observation. 

- In the case of superluminal sources, normalisation of the expected sig- 
nal follows from the most restrictive upper limit on the neutrino sig- 
nal from extraterrestrial sources given by the AMANDA experiment, see 
f lAchterberg et all . l20~07h 



9 dN u 



GeV 



s sr cm" 



(17) 



With an average E~ 2 spectrum for both neutrinos and protons, the spectra 
are connected by assuming that the expected neutrino energy fluence is a 
fraction q of the proton spectrum 



dN u 
dK 



E v dE v = q 



dNp F dF 



with q = 1/40, since only 20% of the proton flux goes into pion production 
via the delta resonance, 1/2 of the remaining flux goes i nto the c harge d 
pion component of which 1/4 goes into neutrinos, see e.g. iBeckerl (120081 ). 
In the case of subluminal sources, using neutrino flux limits leads to 
an excess above the observed spectrum of charged CRs, since the lim- 
its are not stringent enough yet. Instead, the measured part above the 
'ankle' of the CR spectrum is used for an estimate of the contribution 
fro m subluminal sources. The CR en ergy flux above the ankle is given 
by flWaxman and Bahcalll . ll997L Il999h 



] E (E m in = 3 • 10 18 eV) := J 



dK 



v - EL dE n « 10~ 7 GeV cm" - s - sr 



-2 „-l 



3io 18 eV 



dE p v 
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It is expected that this contribution comes from a combined signal from 
AGN and GRBs. In the following it is assumed that the fraction of UHE- 
CRs coming from AGN, contributes a fraction < x < 1. Therefore, the 
fraction of UHECRs from GRBs is (1 - x). 

Thus, the total spectrum as observed at Earth is given as 

^ = A p J (x • ^{E P {z)) + (!-*)■ d -^(E p (z))) ■ g{z) dz . (20) 

•^min 



The minimum redshift is set to z — 0.0018 as the distance of the closest AGN, 
Centaurus A. The maximum redshift is taken to be z max = 7. As the main 
contribution comes from redshifts of z ~ 1 — 2 due to the high number of 
sources at these redshifts, the exact values of the integration limits are not 
crucial. 



3.2.3 Comparison with the observed cosmic ray spectrum 

The diffuse spectrum as measured at Earth is shown in Figure 3. Data points 
represent measurements from a selection of experiments. Our calculated spec- 
tra from superluminal and subluminal shocks are displayed as the dashed and 
solid lines. 

It appears from Figure 3, that in the superluminal shock case, the only possible 
contribution to the measured CR spectrum is around the knee. It is expected, 
however, that the effective flux is actually even lower, because the normalisa- 
tion is based upon the assumption that the contribution cannot be more than 
the current neutrino flux limits permit. Therefore, the calculated flux can be 
considered as an absolute upper limit. Here, the fraction of AGN protons has 
been chosen to be x=0.5, assuming that 50% of the signal is produced by AGN 
and 50% by GRB 

The subluminal shock case has been investigated for different scenarios. The 
upper line represents a spectrum that would be produced by AGN only (x = 
1). The lower line represents a pure GRB spectrum (x = 0). The flux is 
too low to explain the observed component above the ankle if a significant 
contribution comes from GRBs because of the flatness of the GRB spectra. 
It seems that the flat spectra do not fit the present observations. Thus, AGN 
are the favoured sources for the production of UHECRs. Figure 4 shows the 
CR spectrum at the highest energies, multiplied by E 2J in order to have a 
clearer view on the features of the spectrum. The upper and lower solid lines 
represent the same predictions for shocks as in Fig. 3. The middle line assumes 
that 50% is made up by AGN and the remaining 50% comes from GRBs 
The normalisation fits the HiRes data. There is a discrepancy between the 
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Fig. 3. The maximum predicted diffuse flux from GRBs and AGN with superlumi- 
nal shock fronts (dashed line) and subluminal shocks (solid lines). For subluminal 
sources, the upper line is a pure AGN-produced spectrum, the lower line repre- 
sents a pure GRB spectrum. The flux is compared to the measur ed CR spectrum 



Data points are taken from the different experiments: Auger - lYamamoto et al 



(l2007h: HiRes - iThe High Resolution Fly's Eve Collaboration! |200^); AGASA 



Yoshida et al. ( 1995); Yakutsk - Krasilnikov et al. ( 1985 ); Haverah Park - Aye et all 



{2001|); 
Akeno 



HEGRA 



Aharonian et al] (Il999h: C ASA-MIA -iGlasmacher et al.1 (Il999h: 



KrasilnikovelaLl (119831: Tibet - lOzawa et al.1 (120031): Tien Shan 



Antonov et all dl995l\; MSU -iKhristiansen et al.1 dl994h;~J ACEE -lAsakimori et~aD 



Antoni et al 



(|200.5l l In 



dl995j); Proton-Sat - iGrigorov et al.1 (jl975h : KASCADE 
the case of superluminal sources, 50% is assumed to come each from GRB and from 
AGN. 
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Fig. 4. Spectrum of UHECRs m ultiplied by E 2 ' 7 . Data points from Auger 
(lYamamoto et al.1 . 120071 ) and HiRes (jThe High Resolution Fly's Eye Collaboration! . 



20021 ). The solid lines represent the same predictions as presented in Fig. 3. Auger 
data have been renormalised at 10 10 3 GeV to HiRes data and the calculated spec- 
tra have also been normalised to HiRes data. The data can be described well by a 
pure AGN spectrum (blue line) within experimental uncertainties and test particle 
acceleration accuracy. The red line is a mixture of 50% GRB contribution and 50% 
AGN contents, the black line is a pure GRB spectrum. Both do not fit the data. 



normalisation of HiRes and Auger data, which is not entirely understood yet, 
but is probably due to systematic errors in the energy and flux determination 
of the experiments. Therefore, the Auger data are renormalised at 10 10 ' 3 GeV 
to match the normalisation of the HiRes data. With an assumed significant 
contribution from GRBs, it seems difficult to explain the observed spectrum. 
Within the uncertainties of the experimental data, our predicted spectra fit 
the measured spectra between 10 9,8 GeV and 10 10 ' 5 GeV if pure AGN spectra 
are assumed. The flux is slightly too low for the first data points above the 
ankle, between 10 9,5 GeV and 10 9 ' 8 GeV. The steepening of the spectrum at 
lower energies cannot be due to some contribution from sources with a high 
field inclination to the normal because we have shown that such superluminal 
shocks cannot produce energies above about 10 10 ' 5 GeV. It is possible that 
the steepening arises from the addition of subluminal sources with varying 
maximum energy, determined by varying maximum field strengths, so only 
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a frac tion of the sources reach the highest energies of 10 21 eV. lAhlers et al. 



( 120051 ). fitting HiRes data to a cosmological source distribution similar to ours 
in Section 3.2.2, find some steepening around 10 9 GeV due to details of the 
photo-production propagation function. While these effects suffice to explain 
the small difference for the lowest energies above the ankle concerning AGN 
spectra, they cannot explain the large differences between the data and the 
GRB spectral predictions. 



4 Implications for high energy neutrino and photon astronomy 



Fig. 3 has an interesting implication for neutrino and TeV-photon astronomy. 
Neutrinos and TeV photons are pr oduced in proton-photon or proton-nucleon 
interactions, see e.g. iBeckerl (120081 ) for a review, 



pp 



A- 



p7T° , fraction 2/3 
nn + , fraction 1/3 



H + It 7T° 



(21) 
(22) 



The decay of the 7r° leads to high energy photon emission, and n^— particles 
produce neutrinos. The photon signal at TeV energies is not unique, since 
leptonic processes like Inverse Compton scattering contribute at the same 
energies. Therefore, the best, unambiguous way of identifying the hadronic 
interactions are neutrino observations. 



Models for neutr i no em i ssion in AGN a re present in e .g. iMannheiml ( 119951 ) ; 
Mannheim et all (l200lh : ISteckerl tj2005h : iBecker et aD (j2005|), where it is as- 
sumed that protons accelerated in AGN jets can interact with different photon 
fields to produce neutrinos. 

Proton-photon interactions in the prompt phase of GR Bs can lead to neutrino 



and TeV-ph oton production. In the first approach by IWaxman and Bahcall 



( 119971 . Il999l ). it is assumed that GRBs are the sources of UHECRs to cal- 
culate the neutrino spectrum. According to our results in the previous sec- 
tions, GRBs are unlikely to be the sources of UHECRs and the model does 
not hold anymore. However, further developments of the model normalize the 



Guetta et al. 


( 


2004): 


Becker et al. ( 


2006 


); 



e.g. 

Here, the protons do not need to be accelerated to the highest energies: The 
photon field of the prompt emission of GRBs has characteristic energies of 
around ~ 100 keV— 1 M eV, and can reach energies up to > 100 MeV, see 



e.g. (ISchneid et all Il992l : iGonzalez et all . 120031 : iHurlev et all ll99J) . The pro 
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ton energy necessary to produce a Delta resonance, and with this TeV-photons 
and neutrinos, is given as 



E n > 



m, 



mz 



(23) 



Here, tua is the mass of the Delta resonance, m p is the proton's mass, z is the 
redshift of the source and E 1 is the characteristic photon energy. The proton 
energy required for the process is therefore given as 




Thus, with redshifts typically of the order of z = 1 and boost factors of around 
T = 300, the proton energy sufficient for neutrino and TeV-photon production 
is as low as 

£ p ~10 6 GeV. (25) 



This opens the possibility of neutrino production from sources which are not 
observable in charged UHECRs. While the lower energy spectra of extragalac- 
tic, charged cosmic rays cannot be observed due to the high galactic back- 
ground, neutrinos may serve to investigate those sources further. 

Using the same reasoning as above, we can conclude that charged UHECRs 
from superluminal spectra may not be observed in charged cosmic rays, but 
are good candidates for the production of high energy neutrinos and photons 
from AGN or GRBs. Those shocks may produce low energy spectra of much 
higher intensity than subluminal shocks, which are simply hidden due to the 
galactic background. Neutrinos and photons, on the other hand, point back to 
the original source, and may be identified. We can use the maximum energy 
of the proton spectra to calculate the energy of the photon spectra which are 
necessary to produce high energy photons and neutrinos by using Equ. (23), 




For mildly relativistic shocks of V ~ 10 as they occur in AGN, superlumi- 
nal proton spectra reach up to E p ~ 100 GeV. This requires photon en- 
ergies of E 1 > 10 GeV. These high energies can be produced by Inverse 
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Compton scattering of synchrotron or external photons with the accelerated 
electrons. A catalog of AGN with photon emi ssion above 100 MeV w as al- 



ready presented by the EGRET experiment (ISreekumar et al.l . Il998l ). and 



more sources are likely to be id entified when GLAST is launched this summer 
( IGehrels and Michelsonl . Il999l ). 



Highly relativistic shocks with T > 100 reach up to proton energies of 10 5 GeV. 
This requires photon energies of > 100 MeV. This is about 2 orders of mag- 
nitude above the characteristic energy, but the photon spectrum is likely to 
extend to energies aboye 100 MeV as already observed for more than 30 GR Bs, 



see e.g. f lSchneid et all Il992l ; iGonzalez et al.l . 120031 : iHurlev et all . 1199J). If 



mildly relativistic internal GRB shocks are to contribute to neutrino produc- 
tion, photon energies in excess of 1 GeV are required. Such sources would also 
need to be insignificant producers of UHECRs, to satisfy the discussion of 
Section 3.2.2. 



In conclusion, low energy proton spectra from extragalactic sources cannot 
be observed directly due to the high galactic background of cosmic rays, but 
they may be detected indirectly by means of high energy neutrino and photon 
spectra. 



5 Summary & conclusions 



In this work we have presented Monte Carlo simulation studies of the ac- 
celeration of test particles in relativistic, subluminal and superluminal shock 
environments. The source candidates discussed were AGN jets with mildly- 
relativistic shocks of boost factors of T pa 10 — 30 and GRBs regions with 
highly-relativistic shocks, 100 < T < 1000. The resulting particle spectra were 
used to calculate a contribution to the diffuse CR spectrum. 

Particle spectra have been obtained with varying the shock boost factor T and 
shock obliquity, i.e. the inclination angle between the shock normal and the 
magnetic field, ip. Only subluminal shocks are efficient enough to accelerate 
particles up to 10 12 GeV, while superluminal shocks are effective up to ~ 
10 5 GeV. Flat spectra are found for very high subluminal shock boost factors, 
but for superluminal shocks the spectral indices stay roughly constant between 
values of 2.0 to 2.3 in the limited region of efficient acceleration, before a 
cutoff sets in. For the subluminal shock cases, the spectra for mildly-relativistic 
shocks have spectral indices around 2.0 — 2.2. Highly-relativistic shocks have 
spectra as flat as E p ~ - 7 -Ep' 1 - 5 at energies between 10 9 ' 5 GeV and 10 10 - 5 GeV. 
There is no universal spectral form, rather a variety of spectral shapes with a 
noticable plateau-like structure developing at higher T values. This structure is 
very probably related to the number of scattering cycles undergone by particles 
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a particular energy. 



ur results can be summarised as follows: 



Subluminal shock studies with the pitch scattering angle determined to lie 
in the range 1/T < 56 < 10/T, approximately correspond to a situation 
with a spectrum of scattering waves, P(k) B 2 = 5/4/a/2 • T~ 2 ■ k~ x and 
with the neglect of cross-field diffusion. The resulting spectral slopes were 
roughly independent of inclination angle, though details of the features 
were different. A dependence of the spectral index a p on the shock boost 
factor T was found, leading to spectra of a p ~ 2.0 — 2.1 for mildly- 
relativistic shocks of T ~ 10 — 30, but producing much harder spectra 
(0.7 < a p < 1.5) for highly-relativistic shocks, 100 < V < 1000. 

The above implies that GRB particle spectra arising from relativistic 
shocks with very high boost factors between 100 < V < 1000, have spectra 
flatter than E p ~ 1,5 , much flatter than AGN spectra, ~ E p ~ 2 . Moreover, 
the ab ove fi ndings are supported by the lower T work of lNiemec and Ostrowski 
( 120051 ) and lStecker et all (120071 ). Observational evidence, (jDingusl . Il9951 ). 
regarding irregular and flat spectra from GRBs may be explained by the 
spectra we present. This work is also consistent with the general obser- 
vations of the electron spectra that may be injected from the terminal 
hotspots to the lobes of the powerful FR-II radio galaxies which are not of 



a sing le and universal power-law form, as shown in detail in lRudnick et al. 
(119941 ) . iMachalski et al.l (l2007h . etc. 

2) Superluminal shocks are only efficient in accelerating CRs up to E p ~ 
10 5 GeV, resulting in spectral indices of a p ~ 2.0 — 2.3. On the other 
hand, subluminal shocks are more efficient and able to accelerate CRs up 
to E p ~ 10 12 GeV, factors of lO 9 ^ 11 above the particle injection energy. 

3) We discussed the possible contributions of AGN and GRBs to the UHECR 
flux. For superluminal sources, such contributions can be excluded using 
current neutrino flux limits to normalise the spectrum. In the case of 
subluminal sources, the spectrum is normalised to the CR flux above the 
ankle, E min = 10 9 - 5 GeV. Using only AGN (r = 10), the spectrum fits 
the data within experimental uncertainties. With a significant contribu- 
tion from the very high relativistic shocks in GRBs (100 < T < 1000), 
however, the total spectrum is too flat and it is difficult to explain the 
lower part of the spectrum around E p ~ 10 9 5 GeV. Even if UHECRs are 
accelerated in either external or internal shocks of GRBs, it is necessary 
to account for all the energy in accelerated particles, down to the injec- 
tion energy. The total relativistic plasma output available from GRBs is 
only marginally sufficient to account for the total energy required in the 
extra-galactic CR spectrum. 

4) Recent Auger results indicate a correlation between CRs at the highest 

energ ies (E p > 5.7- 10 10 GeV) and the distribution of AGN (IThe Pierre Auger Collaboration 



20071 ). This is a first evidence that the cosmic ray flux above the GZK 
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cutoff originates from AGN predominantly in the supergalactic plane. 
The question of the origin of CRs below the GZK cutoff is not answered 
by this observation, but it is likely that more distant AGN contribute 
significantly to the flux, as AGN in the supergalactic plane make up the 
flux above the GZK cutoff. Moreover the output of X-ray active AGN 
within 60 Mpc provides an energy density exceeding the local estimated 
total extra- galactic CR energy density by a factor 10 4 . It therefore seems 
reasonable to believe that the results from Auger can be explained by 
subluminal relativistic shock acceleration in AGN. It is now important 
to further develop the simulation results by including particle interac- 
tions, more detailed modelling of particle propagation in the inter-galactic 
medium and by investigating the source distribution of AGN in relation 
to observation in order to resolve the question of which AGN are the 
main sources of UHECRs. 
(5) Extragalactic, superluminal shocks are good candidates for the produc- 
tion of high energy neutrinos and photons. The energy density at low 
energies may be quite high compared to the observed flux of UHECRs, 
but hidden by galactic cosmic rays. Neutrino- and photon fluxes may, 
on the other hand, be identified and are probably the only possibility to 
observe extragalactic, superluminal shocks. 
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A The simulation 



The use of a Monte Carlo technique to solve Equ. (5) is dependent on the 
assumption that the collisions represent scattering in pitch angle and that the 
scattering is elastic in the fluid frame where there is no residual electric field. 
Since we assume that the Alfven waves have lower speed than the plasma flow 
itself, the scattering is elastic in the fluid frame. A phase averaged distribution 
function is appropriate to the diffusion approximation we employ which uses 
many small angle scatters. 



We begin the simulation by injecting 10 5 particles far upstream and of a 
weight w p equal to 1.0. A splittin g technique is used similar to the one used 
in the Monte Carlo simulations of iMeli and Quenbv fl2003al ]bh. so that when 
an energy level is reached such that only a few accelerated CRs remain, each 
particle is replaced by a number of N particles of statistical weight 1/N, so 
as to keep a roughly constant number of CRs followed. First order Fermi 
(diffusive) acceleration is then simulated by following the particles' guidance 
centres and allowing for numerous pitch angle scatterings in interaction with 
the assumed magnetised media, while at each shock crossing the particles 
gain an amount of energy determined by a transformation of reference frame. 
The particles are assumed to be relativistic with an initial injection energy 
of 7 ~ (T + 100) when they are entered in the model upstream directed 
towards the shock. As a justification f or a t est particle approach, we note that 
Bell ( 1978al lbl) and Jones and Ellison ( 1991 ) have shown that 'thin' sub-shocks 



appear even in the non-linear regime, so at some energy above the plasma T 
value, the accelerated particles may be dynamically unimportant while they 
re-cross the discontinuity. Another way of arriving at the test-particle regime 
is to inject particles well above the plasma particle energy when they are 
dynamically unimportant and thus require the seed particles to have already 
been pre-accelerated. 



The basic coordinate system employed to describe a shock is a Cartesian 
system (x, y, z), where the shock plane lies on the (y, z) plane. The refer- 
ence frames used during the simulations are the upstream and downstream 
fluid frames, the normal shock frame (NSH) and the de Hoffmann- Teller (HT) 
frame, see Figures A.l and A. 2. 

For the oblique shock cases studied here, provided the field directions encoun- 
tered are reasonably isotropic in the shock frame, we know that tan^i = 
Yi 1 tanip^sH ~ IT 1 ~ ipi where '1' and 'NSH' refer to the upstream and 
normal shock frames respectively. The concentration of field vectors close to 
the x-axis in the upstream fluid frame allows for a reasonable probability of 
finding a de Hoffmann- Teller frame with a boost along the negative y-axis 
less than c. Making this boost then yields an upstream HT frame inclination, 
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Fig. A.l. The coordinate system of a shock as seen in the so called (normal) shock 
frame. 



Fig. A. 2. The coordinate system of a shock as seen in the so called de Hoffmann 
and Teller frame. 

tanif)HT,i — Tht,\ tan^i. While all particles are allowed to cross from down- 
stream to upstream, only particles with a critical HT frame pitch angle, 9 C , 
given by 




-x 



x=0 



+x 



v =0 

V SH 
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are allowed to cross down to upstream and conservation of the first adiabatic 
invariant is used to determine the new, downstream pitch angle. A compres- 
sion ratio of 3 is used al though some MHD conditions favour a value of 4. 



Meli and Quenbyl (j2003bl ) do not find a consid erable differenc e in sim ulation 



results between these two cases. The results of iNewman et all ( 119921 ) suggest 
that it is legitimate to use Equ. A.l as a reasonable approximation. These 
authors checked the preservation of the first adiabatic invariant for the worst 
case, near perpendicular shock, employing trajectory integration with a re- 
alistic scattering field right up to the shock interface. At the critical angle 
for a compression ratio of three shock, adiabatic invariant deviations were 
typically confined within ten percent while large effects tended to occur only 
towards 90 degree pitch angle. Particles are assigned a random phase so that a 
3-dimensional transformation of momentum vectors can be achieved between 
the fluid and HT frames. Away from the shock, the guiding centre approxi- 
mation is used so that a test particle moving a distance, d, along a field line 
at xjj to the shock normal, in the plasma frame has a probability of collision 
within d given by P(d) = 1 — exp(— d/X) = R, where the random number R 
is < R < 1. Weighting the probability by the current in the field direction 
H yields d = —XfilnR. The pitch angle is measured in the local fluid frame, 
while the value Xi gives the distance of the particles to the shock front, where 
the shock is assumed to be placed at x = 0. Furthermore, Xi is defined in 
the shock rest frame and the model assumes variability in only one spatial 
dimension. Scattering in pitch is applied as described in the main body of this 
text. 



In our simulations continuous Lorentz transformations are performed from 
and into the local plasma frames into or from the shock frame in order to 
check for particle shock crossings. All particles leave the system if they escape 
far downstream at the spatial boundary, r b . The downstream spatial bound- 
ary required can be estimated initially from the solution of the convection- 
diffusion equation in a non-relativistic, large-angle scattering approximation 
in the downstream plasma, which gives the chance of return to the shock as 
exp(— V^rb/xj). In fact, we have performed many runs with different spatial 
boundaries to investigate the effect of the size of the acceleration region on 
the spectrum, so as to find a region where the spectrum is size independent. 
Alternatively, the particles leave the system if they reach a specified maximum 
energy E m&x for computational convenience. 



For the superluminal shock conditions, where the physical picture of the shock 
drift acceleration applies, it is necessary to abandon the guiding centre ap- 
proximation when the trajectories begin to intersect the shock surface. Here, 
we consider a helical trajectory motion of each test-particle of momentum p, 
in the fluid frame, upstream or downstream, where the velocity coordinates 
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(v x ,v y ,v z ) of the particle are calculated in 3-dimensional space as follows 
v Xi = Vi cos 9i cos fa — Vi sin 9i cos fa sin fa , (A. 2) 

v Vi = V i COS @i + V i COS 4>i cos V't (A -3) 

and 

v Zi = -Vi sin ^ sin fa . (A.4) 

Here, 0j is the pitch angle, 0, G (0, 27r) and ^ is the angle between the magnetic 
field and the shock normal in the respective fluid frames (i—1,2 for upstream 
and downstream respectively). 

We follow the trajectory in time, using fa = <p + ut, where t is the time from 
detecting shock presence at xnsh, Vnsh, znsh by using 

dx = x NS H + v Xi 5t , (A. 5) 

dy = Vnsh + v Vi 5t (A. 6) 



and 

dz = z NSH + v z M (A. 7) 



assuming that St = r g /Hc where H > 100 and r g is the Larmor radius. The 
particle's gyrofrequency to is given by the relation, uji = e|-Bj|/7j, Bi is the 
magnetic field, 7, is the particle's boost factor and e is its charge in gaussian 
units. 

We follow the helical trajectory of each particle in time t, in the new frame 
where t is the time from detecting the shock intersection at (x, y, z) until the 
trajectory has performed one gyro period without re-intersecting the shock 
surface. Nevertheless, because of the peculiar properties of the helix we need to 
establish where a particle, starting off in the upstream frame, with a particular 
9 and <fi first encounters the shock. To establish when the shock encounter 
happens, we choose to go back a whole period, Tj = 2-K/uji by reversing signs of 
the helix velocity coordinates and by keep checking throughout the simulation 
to determine if the particle trajectory encounters the shock front, placed at 
x = in the shock rest frame. If the particle encounters the shock then the 
suitable Lorentz transformation to the relevant fluid rest frame is made and 
we continue following the particle helical trajectory until shock intersections 
cease. At this juncture, the guiding centre is followed in the same way as in 
the diffusive acceleration picture of the subluminal shocks. During the helical 
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phase of the numerical integration, the prescription for pitch angle scatter is 
applied as in the general plasma frame motions, in order to more realistically 
simulate a mean plus chaotic field situation where turbulence is clearly present 
close to a shock. 



B The spectra 



The particle spectra resulting from the simulation are presented in Figures B.l, 
B.2 and B.3. In each of the figures, the simulated particle spectra are shown 
(black dots) for three different shock angles each, ip = 23°, 33°, 43°, for nine 
different boost factors, T = 10, 20, 30 in the first row, starting from the left, 
r = 100, 300, 500 in the middle row and T = 700, 900, 1000 in the lower row. 
Each graph shows the logarithm of the proton spectrum <i$ / dE p in arbitrary 
units versus the proton's energy in units of GeV. Note that the spectrum 
can be expressed more generally in terms of the particle's boost factor 7 = 



E p /(m p 



Therefore, the results are also valid for nuclei with higher mass 



(e.g. Fe). In the present investigation protons are considered. 

Figures B.l, B.2 and B.3 show that the spectral shape starts to deviate from 
a power-law as T increases with the onset of plateau formation. This might be 
expected, since the particles are swept away rapidly downstream with a low 
possibility to return upstream for high T factors (around 20% of the particles 
return to the shock after one shock cycle). A small chance to return to the 
shock (except for a relatively small subset of downstream pitch angle particle 
'histories') produces the spectral irregularities and the anisotropy seen for 
these returning particles. We note that in the cases of T = 10 — 30 relatively 
smooth spectra are produced, but spectra become more structured (plateau- 
like) at the more extreme values T — > 1000. In the latter cases the effects of 
individual acceleration cycles are clearly evident. The mechanism of plateau 
develo pment as an acc eleration cycle effect is implicit in figure 6 of IProtheroe 
(|200ll ) and figure 2 of IStecker et al.l (120071 ) and seems to be independent of 
the shock inclination angle for a particular scattering model. The structured 
sp ectra can also be seen in the l ower T large angle s catter model simu l ations 
of buenbv and Lieul fll989h and lEUison et all fll990h while IProtheroe! (l200lh 
shows a similar contrast in behaviour, between large and small angle scattering 
models up to T = 20. 



The initial bump in most of the spectra is due to the monochromatic injection 
of the particles. We approximate the spectra with a power-law (straight lines 
seen in the figures) at higher energies and thus, artificial injection features are 
excluded from the calculations. The maximum energy is chosen to be 10 21 eV 
as discussed in Section 1.2. The plateau-like parts in the spectra at higher 
energies are physical features, especially at high boost factors, as particles 
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continue to be accelerated in a second cycle. The lower energy part of the 
spectrum is dominated by particles undergoing one acceleration shock cycle, 
while the second bump in the spectrum represents particles having completed 
two cycles. A shock cycle is a shock crossing from upstream to downstream 
to upstream. Since some upstream particles can suffer reflection before down- 
stream transmission, this effect increases the statistical energy gain of the 
particles in their overall encounter with the shock surface. The aim of fitting 
a single power-law to the computed spectral points is primarily to examine 
the variation of primary spectra with the T factor of the shock. The spectral 
points themselves are difficult to compare because the structure becomes more 
complex with the increasing boost factor. As discussed in the main text, the 
relevant way to compare the particle spectra to the data is to fit the energy 
range observed in cosmic rays. The straight lines in Figures B.l, B.2 and B.3 
together with the values written in the lower left corner of each graph indicate 
the single power-law for comparison chosen individually for each case. 

Table B.l shows the spectral index dependence of the boost factor for the three 
shock inclination angles ip = 23°, 33°, 43°. For each of the three angles, the 
spectra become harder with the increasing T factor of the shock. Due to the 
structure in the spectra, these values represent simple first order approxima- 
tions which are useful for comparison of the simulation results with the data. 
CR data include large statistical and systematic errors, which would make 
it difficult to distinguish features attributable to the acceleration mechanism 
other than single or broken power-laws. 

The spectral flattening with higher T, implies that GRB particle spectra, 
arising from relativistic shocks with very high boost factors between 100 < 
T < 1000, have spectral indices ranging between a p ~ 2.1 — 1.5. These spectra 
tend to be somewhat steeper than some particular entries in Table B.l, but the 
tendency of a flattening with the boost facto r is always present. T he relativistic 



flatten ing: effect is consistent with studies of [BaringJ (119991 . 120041 ) ; IStecker et al 



( 120071 ). On the other hand, for shocks with T ~ 10 — 30 occurring in AGN, 
spectral indices have values between 2.0 < a p < 2.3. 

In addition, as one sees from Table B.l, the single power- law comparison gives 
comparable values for different angles ip. Test simulation runs we performed, 
using a series of values of ip, follow the same trend as long as the scatter- 
ing model i s fixe d. The present results deviate f rom earlier investigations by 
Kirk et al. ( 2000l ); Ostrowski and Bednar2 ( 2002 ). which show a saturation of 



the hardening of the spectra with the boost factor, a p —> 2.33 for r — > oo. 
Nevertheless, the latter studies concern the special cases of extremely small 
values for pitch angle scattering, for particle acceleration in relativistic shocks. 
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Fig. B.l. Subluminal spectra for ip = 23° and different T: T = 10, 20, 30 is displayed 
in the first row, in the middle, T = 100, 300, 500 is shown and T = 700, 900, 1000 is 
the bottom row. The black dots in each graph represent the simulation result. The 
straight line shows the single power-law for comparison. The spectral behaviour is 
indicated in the lower left corner of each graph. 
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Fig. B.2. Subluminal spectra for ip = 33° and different T, as in Fig. B.l. 
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Fig. B.3. Subluminal spectra for ip = 43° and different T, as in Fig. B. 



r a p {ip = 23°) 


ap(</> = 33°) 


a p (V> = 43°) 


10 2.1 


2.1 


2.3 


20 2.0 


2.0 


2.3 


30 2.1 


2.0 


2.2 


100 1.8 


1.8 


2.2 


300 2.0 


1.8 


2.0 


500 1.9 


1.7 


1.6 


700 1.8 


1.4 


1.7 


900 1.5 


1.0 


1.3 


1000 1.2 


1.2 


1.5 



Table B.l 

Spectral indices for a single power-law comparison for subluminal shocks of different 
boost factors and three inclination angles. 
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